function [Ak, phik, FCCoverage, ak, bk] = funcCellContourFCManual(contourX, contourY, noOfFCs)
% DFT on contour (polygon)
% 
% contourX & contourY - x & y coordinates of open polygon as column vectors
% noOfFCs - number of Fourier coefficients to be calculated including a_0
% 
% ak - coefficients a_k
% bk - coefficients b_k
% Ak - overall amplitude: Ak = sqrt(a_k^2 + b_k^2)
% phik - corresponding phase: phik = atan2(a_k, b_k)
% FCCoverage - quality measure stemming from Parseval's theorem

%% Contour
geom = polygeom(contourX, contourY);
centerX = geom(2); centerY = geom(3);
contourX = [contourX; contourX(1)];
contourY = [contourY; contourY(1)];
% origContour = [contourX contourY];

%% Interpolation (linear in xy)
t = linspace(0,1,51); t(end) = [];
tmp = repmat(contourX(1:end-1),1,length(t))+(contourX(2:end)-contourX(1:end-1))*t; tmp = tmp'; tmp = tmp(:);
interpX = tmp;
tmp = repmat(contourY(1:end-1),1,length(t))+(contourY(2:end)-contourY(1:end-1))*t; tmp = tmp'; tmp = tmp(:);
interpY = tmp;

%% Transformation in polar coordinates (including sorting)
% origR = sqrt((contourX - centerX).^2 + (contourY - centerY).^2);
% origPhi = atan2(contourY - centerY, contourX - centerX); origPhi = mod(origPhi, 2*pi);
% a = sortrows([origPhi(1:end-1) origR(1:end-1)]);
% origPhi = a(:,1); origR = a(:,2);

interpR = sqrt((interpX - centerX).^2 + (interpY - centerY).^2);
interpPhi = atan2(interpY - centerY, interpX - centerX); interpPhi = mod(interpPhi, 2*pi);
a = sortrows([interpPhi interpR]);
interpPhi = a(:,1); interpR = a(:,2);

%% Fourier Transformation
interpdPhi = mod(diff([interpPhi; interpPhi(1)]),2*pi);
interpdPhi = 0.5*(interpdPhi + [interpdPhi(end); interpdPhi(1:end-1)]);

A0 = 1/pi * sum(interpR .* interpdPhi);
ak = 1/pi * sum(repmat(interpR,1,noOfFCs) .* cos(repmat(0:noOfFCs-1,length(interpR),1) .* repmat(interpPhi,1,noOfFCs)) .* repmat(interpdPhi,1,noOfFCs) ,1);
bk = 1/pi * sum(repmat(interpR,1,noOfFCs) .* sin(repmat(0:noOfFCs-1,length(interpR),1) .* repmat(interpPhi,1,noOfFCs)) .* repmat(interpdPhi,1,noOfFCs) ,1);

Ak = sqrt(ak.^2 + bk.^2);
phik = atan2(ak, bk);

Ak(1) = A0/2;
ak(1) = 0.5*ak(1);
%% Fourier reconstruction
% rRec = A0/2 + sum(repmat(Ak(2:end),length(phi),1) .* sin(repmat(1:noOfFCs-1,length(phi),1) .* repmat(phi,1,noOfFCs-1) + repmat(phik(2:end),length(phi),1)),2);
% xRec = rRec .* cos(phi);
% yRec = rRec .* sin(phi);

%% quality measure
parsevalAll = 1/pi * sum(interpR.^2 .* interpdPhi ,1);
parsevalNFC = A0^2/2 + sum(Ak(2:end).^2);
FCCoverage = parsevalNFC / parsevalAll;

end